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In this paper, we explore the interior dynamics of neutral and charged black holes in f{R) gravity. 
We transform f{R) gravity from the Jordan frame into the Einstein frame and simulate scalar 
collapses in flat, Schwarzschild, and Reissner-Nordstrom geometries. In simulating scalar collapses 
in Schwarzschild and Reissner-Nordstrom geometries, Kruskal and Kruskal-like coordinates are used, 
respectively, with the presence of f' and a physical scalar field being taken into account. The 
dynamics in the vicinities of the central singularity of a Schwarzschild black hole and of the inner 
horizon of a Reissner-Nordstrom black hole is examined. Approximate analytic solutions for different 
types of collapses are partially obtained. The scalar degree of freedom (j>, transformed from /', plays 
a similar role as a physical scalar field in general relativity. Regarding the physical scalar held in 
f{R) case, when dcp/dt is negative (positive), the physical scalar field is suppressed (magnified) by 4>, 
where t is the coordinate time. For dark energy f{R) gravity, inside black holes, gravity can easily 
push /' to 1. Consequently, the Ricci scalar R becomes singular, and the numerical simulation 
breaks down. This singularity problem can be avoided by adding an R^ term to the original f{R) 
function, in which case an infinite Ricci scalar is pushed to regions where f' is also infinite. On the 
other hand, in collapse for this combined model, a black hole, including a central singularity, can 
be formed. Moreover, under certain initial conditions, f' and R can be pushed to infinity as the 
central singularity is approached. Therefore, the classical singularity problem, which is present in 
general relativity, remains in collapse for this combined model. 


I. INTRODUCTION 

The internal structure of black holes and spacetime sin¬ 
gularities are key topics in gravitation and cosmology dl¬ 
ls], and are great platforms to explore the connection 
between classical and quantum physics. It is widely be¬ 
lieved that rotating black holes exist in reality. According 
to Price’s theorem, for a collapsing star, the gravitational 
radiation carries away all the initial features of the star’s 
gravitational field, except the mass, charge, and angular 
momentum parameters [Bj. As a further step, it is nat¬ 
ural to ask what the final state of the internal collapses 
might be. 

Ever since the foundation of general relativity, people 
have been trying to go beyond it. This endeavor arises 
from unifying gravitation and quantum mechanics, and 
addressing some cosmological problems, including the 
singularity problem in the early Universe and the dark 
energy problem in the late Universe. Various modihed 
gravity theories have been explored, including scalar- 
tensor theory, high-dimensional theory, and f{R) grav¬ 
ity, etc. For a review of modified gravity theories, see 
Ref. [7|. For reviews of f{R) theory, see Refs. [HUTT]. 

Static and spherically symmetric black hole solutions 
in f{R) gravity were explored in Ref. [Tl|. Charged 
Born-Infeld black holes for f{R) theories were stud¬ 
ied in Ref. [13]. Instabilities and (anti)-evaporation of 
Schwarzschild-de Sitter and Reissner-Nordstrom black 
holes in modified gravity were discussed in Refs. [HHis]. 
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Gravitational collapses in some modified gravity theor¬ 
ies have been studied numerically. Spherical collapse of 
a neutral scalar field in a given spherical, charged black 
hole in Brans-Dicke theory was investigated in Ref. [Hi- 
Spherical collapses of a charged scalar field in dilaton 
gravity and f{R) gravity were explored in Refs. |18j 
and m, respectively. Spherical scalar collapse in f{R) 
gravity was simulated in Ref. [20]. Asymptotic analysis 
was implemented in the vicinity of the singularity of a 
formed black hole. 


A. Mass inflation 

In the vicinity of the central singularity inside a 
Schwarzschild black hole, the tidal force diverges, and 
maximal globally hyperbolic region defined by initial 
data is inextendible. However, inside charged (Reissner- 
Nordstrom) and rotating (Kerr) black holes, the central 
singularity is timelike. The globally hyperbolic region 
is up to the Cauchy horizon, and the spacetime is ex¬ 
tendible beyond this horizon to a larger manifold. The 
Reissner-Nordstrom inner (Cauchy) horizon is a surface 
of infinite blueshift, which in turn may cause the inner 
horizon unstable m- Furthermore, the strong cosmic 
censorship conjecture was proposed, which states that 
for generic asymptotically flat initial data, the maximal 
Cauchy development is future inextendible. For math¬ 
ematical explorations of the internal structures of charge 
black holes, see Refs. [HJI^S]. For reviews on the Cauchy 
problem in general relativity and strong cosmic censor¬ 
ship, see Refs. [SUES], respectively. 

The backreaction of the radiative tail from a grav¬ 
itational collapse on the inner horizon of a Reissner- 
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Nordstrom black hole was investigated by Poisson and 
Israel [Ml EH- It was shown that due to the divergence 
of the tail’s energy density occurring on the inner hori¬ 
zon, the effective internal gravitational-mass parameter 
becomes unbounded. This phenomenon is usually called 
mass inflation. These arguments were extended to the 
rotating black hole case in Ref. [55] ■ 

In Refs. [Ml EZ]) approximate analytic expressions 
were obtained by considering a simplified model in which 
the perturbations were modeled by cross-flowing radial 
streams of infalling and outgoing lightlike particles. To 
get more information, some numerical simulations in 
more realistic models have been performed. The dynam¬ 
ics of a spherical, charged black hole perturbed nonlin- 
early by a self-gravitating massless scalar field was nu¬ 
merically studied in Refs. [55H51] . Under the influence 
of the scalar field, the inner horizon of a charged black 
hole contracts to zero volume, and the center becomes 
a spacelike singularity. The mass inflation phenomenon 
was observed. In Refs. [MIES], with regular initial data, 
spherical collapse of a charged scalar field was simulated. 
An apparent horizon was formed. A null, weak mass- 
inflation singularity along the Cauchy horizon and a fi¬ 
nal, spacelike, central singularity were obtained. Spher¬ 
ical collapses in Brans-Dicke theory, dilaton gravity, and 
f{R) gravity were investigated in Refs. [iTUIfl] . Mass 
inflation phenomena were also reported. 

It is important to connect approximate analytic can¬ 
didate expressions with numerical results. In Refs. [Ml 
[88] . the features of the Cauchy horizon singularity in 
charge scattering were studied. Analytic and numerical 
results were compared at some steps. 

In this paper, we use the following notations: 

(i) Neutral collapse: neutral scalar collapse toward a 
black hole formation. 

(ii) Neutral scattering: neutral scalar collapse in a 
(neutral) Schwarzschild geometry. 

(iii) Charge scattering: neutral scalar collapse in a 
(charged) Reissner-Nordstrom geometry. In this 
process, the scalar field is scattered by the inner 
horizon of a Reissner-Nordstrom black hole. 


B. New results 

In this paper, we explore neutral scalar collapses in 
flat, Schwarzschild, and Reissner-Nordstrom geometries 
in /(i?) gravity, taking the well-known Hu-Sawicki model 
as an example [39] ■ We seek approximate analytic solu¬ 
tions. A generalized Misner-Sharp energy in f{R) grav¬ 
ity in the Jordan frame was defined in Ref. m- In this 
paper, for ease of operation, we mainly work in the Ein¬ 
stein frame and compute the Misner-Sharp mass function 
of the general relativity version, instead. Moreover, we 
will investigate a dark energy f{R) singularity problem. 


We explore scalar collapses in both general relativity 
and f{R) gravity. For convenience, we transform the dy¬ 
namical system in f{R) gravity from the Jordan frame 
into the Einstein frame. In the new system, we equi¬ 
valently work in Einstein gravity to which a scalar de¬ 
gree of freedom </>(= a/3/2 In /'/VSttG) and a physical 
scalar field 'i/f are coupled. Basically, cj) plays a similar 
role as what a physical scalar field ip does in Einstein 
gravity. While in f{R) gravity, the physical scalar field 
Ip is suppressed (magnified) when dcp/dt is negative (pos¬ 
itive), where t is the coordinate time. For simplicity, the 
results in general relativity are presented in a separate 
paper mi, and we focus on collapse in /(i?) gravity in 
this paper. 

According to the strength of the scalar field, charge 
scattering can be classified into five types as follows: 

(i) Type I: spacelike scattering. When the scalar field 
is very strong, the inner horizon can contract to 
zero volume rapidly, and the central singularity be¬ 
comes spacelike. The dynamics near the spacelike 
singularity is similar to that in neutral collapse. 

(ii) Type II: null scattering. When the scalar field is in¬ 
termediate, the inner horizon can contract to a place 
close to the center or reach the center. For each 
variable (the metric elements and physical scalar 
field), the spatial and temporal derivatives are al¬ 
most equal. In the case of the center being reached, 
the central singularity is null. This type has two 
stages: early/slow and late/fast. In the early stage, 
the inner horizon contracts slowly, and the scalar 
field also varies slowly. In the late stage, the in¬ 
ner horizon contracts quickly, and the dynamics is 
similar to that in the spacelike scattering case. 

(iii) Type III: critical scattering. This case is on the edge 
between the above two cases. The central singular¬ 
ity becomes null. 

(iv) Type IV: weak scattering. When the scalar field 
is very weak, the inner horizon contracts but not 
much. Then the central singularity remains time¬ 
like. 

(v) Type V: tiny scattering. When the scalar field is 
very tiny, the influence of the scalar field on the 
internal geometry is negligible. 

In this paper, we will explore the dynamics of Types I, 
II, and IV, and obtain approximate analytic solutions for 
the first two. 

By comparing the dynamics in a Reissner-Nordstrom 
geometry and charge scattering, we investigate the causes 
of mass inflation and seek further approximate ana¬ 
lytic solutions with the following improvements. Usu¬ 
ally, double-null coordinates are used in studies of mass 
inflation in spherical symmetry. In the line element of 
double-null coordinates, the two null coordinates u and 
V are present in the form of product dudv. In the equa¬ 
tions of motion, mixed derivatives of u and v are present 
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quite often. In this paper, we use a slightly modified 
line element, in which one coordinate is timelike and the 
rest are spacelike. In this case, in the equations of mo¬ 
tion, spatial and temporal derivatives are usually separ¬ 
ated. This simplifies the numerical formalism and helps 
to obtain approximate analytic solutions. In addition, 
we compare numerical results and approximate analytic 
solutions closely at each step. We compare the dynamics 
for Schwarzschild black holes, Reissner-Nordstrom black 
holes, neutral collapse, and charge scattering. We treat 
the system as a mathematical dynamical system rather 
than a physical one, examining the contributions from all 
the terms in the equations of motion. 

In Ref. [19] where spherical charged scalar collapse in 
f{R) gravity was simulated, a singularity problem was 
reported. When a dark energy f{R) model is used, /' 
can be pushed to 1 easily. Correspondingly, the Ricci 
scalar R becomes singular. This singularity problem can 
be avoided when an R^ model is used instead. However, 
the causes of this singularity problem were not explained. 
In this paper, we will consider a simpler case. Instead of 
simulating the collapse of a charged scalar field, we study 
neutral scalar collapse in a Reissner-Nordstrom geometry. 
The same singularity problem is found. By analyzing 
the contributions from all the terms in the equations of 
motion for (j){=^3/2\nf'/\/8 ttG) with f'=df/dR, we in¬ 
terpret the causes for the singularity problem. Basically, 
near the inner horizon, in the equation of motion for cj), 
the scalar field (f) and the geometry construct a positive 
feedback system. Depending on initial conditions, (j) can 
be accelerated either in positive or in negative directions, 
until singularities are met. In the negative case, (j) can be 
accelerated to negative infinity. Correspondingly, /' goes 
to zero as the central singularity is approached. How¬ 
ever, in the positive case, (j) can be pushed to zero in a 
short time. Correspondingly, /' and the Ricci scalar R 
are pushed to 1 and infinity, respectively. This is the 
cause of the singularity problem. Taking into account 
quantum-gravitational effects at high curvature scale, one 
may obtain an additional R^ term to the Lagrangian for 
gravity [111133]. When this term is added to the f{R) 
function, a singular R is pushed to regions where f' is 
also singular. Therefore, the singularity problem can be 
avoided [42H49] . 

Although the dark energy f{R) singularity problem is 
avoided in the combined model (a combination of a dark 
energy f{R) model and the R^ model), the classical sin¬ 
gularity problem, which is present in general relativity, 
remains in collapse for this model. Under certain initial 
conditions, near the central singularity, d(j)/dt can be pos¬ 
itive. Then the positive feedback system in the equation 
of motion for (j) can push (/) and R to positive infinity. 

This paper is organized as follows. In Sec. [TTJ we build 
the framework for charge scattering, including action for 
charge scattering, the coordinate system, and the f{R) 
model. In Sec. |HI[ we set up the numerical formalism for 
charge scattering. In Secs. |IV|[V] andjVI] scalar collapses 
in flat, Schwarzschild, and Reissner-Nordstrom geomet¬ 


ries will be explored, respectively. In Sec. VH 


we con¬ 


sider weak charge scattering. In Sec. jVIH we discuss 
the causes and avoidance of the singularity problem. In 
Sec. |IX[ the results will be summarized. 

In this paper, we set G = c = 47reo = I. 


II. FRAMEWORK 

In this section, we build the framework for charge scat¬ 
tering in /(i?) gravity, in which a self-gravitating mass¬ 
less scalar field collapses in a Reissner-Nordstrom geo¬ 
metry in f{R) gravity. Compared to general relativity, 
in this process, there is one extra scalar degree of free¬ 
dom f'=df/dR. For convenience, f{R) gravity is trans¬ 
formed from the Jordan frame into the Einstein frame. 
For comparison and verification considerations, we use 
Kruskal-like coordinates, and set up the initial conditions 
by modifying those in a Reissner-Nordstrom geometry 
with a physical scalar field, a scalar degree of freedom /', 
and the potential for /'. The Hu-Sawicki model is used 
as an example. 


A. Action 


The action for charge scattering in f{R) gravity can 
be written as follows: 



'ML 

IGttG 


J- + Lf , 


( 1 ) 


with 




( 2 ) 


Cp 


4 


(3) 


/(i?)/(167rG), and Lp are the Lagrange densities 
for f{R) gravity, a physical scalar field "0, and the elec¬ 
tric field for a Reissner-Nordstrom black hole, respect¬ 
ively. f{R) is a certain function of the Ricci scalar R, 
and G is the Newtonian gravitational constant. F^^u is 
the electromagnetic-field tensor for the electric field of a 
Reissner-Nordstrom black hole. 

The energy-momentum tensor for the massless scalar 
field Ip is 


rTi{lP) 


2 


(4) 


The electric field of a Reissner-Nordstrom black hole can 
be treated as a static electric field of a point charge of 
strength q sitting at the origin r = 0. In the Reissner- 
Nordstrom metric, the only nonvanishing components of 
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are Ftr = —Ftr = —qlr^. The corresponding energy- 
momentum tensor for the electric field is 1501 


= - 


2 

V\9\ K 


4 ^ K ^^9 


T 


STrr^ 


'diag(-l,- 1 , 1 , 1 ). 


(5) 


Although Eq. ([^ is obtained in the Reissner-Nordstrom 
metric, it is valid in any coordinate system, since as seen 
by static observers, the electromagnetic field should be 
purely electric and radial Emin]- 


B. fiR) theory 

The equivalent of the Einstein equation in /(i?) gravity 
reads 

f'Rf.. - - (V^V,. - g^^n) f = (6) 

where f'=df /dR and □ = Vo,V“. The trace of Eq. ([^ 
describes the dynamics of /', 


□r-^^^-^T = o, 

3 3 


(7) 


where T is the trace of the stress-energy tensor De¬ 
fining a new variable x by 


X = 


_ ^ 

di?’ 


and a potential U (y) by 


rrU 2/ - fR 


one can 


rewrite Eq. 0 


as 


Stt ^ 


□X - C7'(x) - y T = 0. 


( 8 ) 


(9) 


( 10 ) 


The field equations for f{R) gravity ^ are somewhat 
different from the more familiar corresponding ones in 
general relativity. Therefore, for convenience, we trans¬ 
form f{R) gravity from the current frame, which is usu¬ 
ally called the Jordan frame, into the Einstein frame, in 
which the formalism can be formally treated as Einstein 
gravity coupled to a scalar field. 

Rescaling y by 


( 11 ) 


k 4 i = \/ 2 Inx. 


one obtains the corresponding action of f{R) gravity in 
the Einstein frame [5] 


Se = 


d*3 




d^xCM 


( 12 ) 


where k = VSttG, 5 ^, 1 . = X ’ = ixR - 

/)/( 2 «;^X^), and a tilde denotes that the quantities are 
in the Einstein frame. The Einstein field equations are 


G^l, — K 


FW I F{M) 

fliy ' ^1/ 


where 


Tjj.t'’ = df,(l)d^(j) - ^g°'l^da(t)df!cj) + V(<f>) 


(M) 


9p{M) _ 
FM) 


(13) 

(14) 

(15) 


Till, is the ordinary energy-momentum tensor for the 
physical matter fields in terms of g^n in the Jordan frame. 
With the expression for the energy-momentum tensor for 
the scalar field ^/> in the Jordan frame, shown in Eq. Q, 
the corresponding expression in the Einstein frame can 
be written as 


= ^ - ^g^vg°‘^da'il3dp'il) 

which gives 


y(b) _ — _ 




X 


X 


(16) 


(17) 


In the Jordan frame, in any coordinate system, the 
energy-momentum tensor for the static electric field of a 
point charge of strength q sitting at the origin r = 0 can 
be expressed as [see Eq. 


7^^^^() = ^-diag(-l,-l,l,l). 


87rrjp4 

Then we have in the Einstein frame, 
T{q)^ = 

ly — y ^ociy 

(9) 


(18) 


g^la JV 

X X 
9^ 

STTX^rjp"^ 


diag(-l,- 1 , 1 , 1 ) 


g' 

87rrEF 


4 •diag(-l,-l,l,l), (19) 


where rjp and Tef are the quantity r in the Jordan and 
Einstein frames, respectively. Since we mainly work in 
the Einstein frame in this paper, we simply use r for 
Tef- We denote the total energy-momentum tensor for 
the source fields as 


2^(total)/i _ 


( 20 ) 


The equations of motion for (f) and ip can be derived 
from the Lagrange equations as 


1 


□((.-E'(<^)- — = 0 , 

V 6 


( 21 ) 
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Figure 1: Initial and boundary conditions for charge scat¬ 
tering. Initial slice is at t = 0. Definition domain for x is 
[—Xb Xb]. u = {t — x)/2 and v = {t + a:)/2. 



X 


Figure 2: Contour lines for r defined by Eq. (281 in a Reissner- 
Nordstrom geometry with m = 1 and q — 0.7. Although the 
exact inner horizon is at regions where uv and are 

infinite, r can be very close to the inner horizon r = r- even 
when uv and take moderate values. 


□■0 - \/ 3 = 0 . 


( 22 ) 


Alternatively, Eqs. (21) and (22) can be obtained from 
the corresponding ones in the Jordan frame. Some details 
are given in the Appendix. 

In the Einstein frame, the potential for (j) can be writ¬ 
ten as 


Vi^) = 


xR- f 

2k^X^ 


Then we have 


V'id)] = — .^ = ~ 

dx dfj) 


(23) 


(24) 


C. Coordinate system 


In the studies of mass inflation, the double-null co¬ 


ordinates described by Eq. (25) are usually used, 
ds^ = —Ae~'^'^dudv + r'^dfl^, 


(25) 


where a and r are functions of the coordinates u and 
V. u and v are outgoing and ingoing characteristics (tra¬ 
jectories of photons), respectively. Eor convenience, in 


this paper, we use a slightly modified form described 
by Eq. (26), obtained by defining u = {t — x)/2 and 
V = {t + x)/2 [5T] . 


ds^ = e {—dt^ + dx'^) + r^dXl^. 


(26) 


This set of coordinates is illustrated in Fig. Similar to 
the Schwarzschild metric, the Reissner-Nordstrom metric 
can be expressed in Kruskal-like coordinates [S2] (also see 
Refs. [121 H?! [Sni US]). So for ease and intuitiveness, we 
set the initial conditions close to those of the Reissner- 
Nordstrom metric in Kruskal-like coordinates, taking into 
account the presence of a physical scalar field ^/>, a scalar 
degree of freedom /', and the potential U{f'). 

In the form of Ref. [53], the Reissner-Nordstrom metric 
in Kruskal-like coordinates in the region of r > r_ can 
be written as 


ds^ = ^ 6 - 2 '=+’' {-dt^+dx^)+r^dn^ 

(27) 

where r±(= mzL — q^) and k±[= {r± — r;p]/(2r^)] 
are the locations and surface gravities for the outer and 
inner horizons of a Reissner-Nordstrom black hole, re¬ 
spectively. r(t,x) is defined implicitly below [53] . 


4uv = f - x'^ = 6 ^''+’’ 




|fc-i 


• (28) 


In this set of coordinates, as implied by Eq. (28), the 
exact inner horizon is at regions where uv and (t — x"^) 
are infinite. However, it is found that, even when uv and 
— x"^) take moderate values, r still can be very close to 
the inner horizon, e.g., r = (1 -|- 10“^°)r_. (See Fig. |^) 
Therefore, at such regions, the interaction between the 
scalar fields and the inner horizon still can be very strong, 
then we can investigate mass inflation numerically. 

This formalism has several advantages as follows: 


(i) In the line element ( [2^ , one coordinate is timelike 
and the rest are spacelike. This is a conventional 
setup. It is more convenient and more intuitive 
to use this set of coordinates. For the set of co¬ 


ordinates described by Eq. (25), in the equations 


of motion, many terms are mixed derivatives of u 
and v\ while for the set of coordinates described by 
Eq. (26), in the equations of motion, spatial and 


temporal derivatives are usually separated. 


(ii) We set initial conditions close to those in the 
Reissner-Nordstrom metric. Consequently, with the 
terms related to the scalar fields being removed, we 
can test our code by comparing the numerical res¬ 
ults to the analytic ones in the Reissner-Nordstrom 
case conveniently. Moreover, by comparing dy¬ 
namics for scalar collapse to that in the Reissner- 
Nordstrom case, we can obtain intuitions on how 
the scalar fields affect the geometry. 
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(iii) The interactions between scalar fields and the geo¬ 
metry are local effects. In Refs. [29l[30], the space 
between the inner and outer horizons are compacti- 
fied into finite space. This overcompactification, at 
least to us, makes it a bit hard to understand the 
dynamics. In the configuration that we choose, the 
space is partially compactified, and the picture of 
charge scattering turns out to be simpler. 


III. NUMERICAL SETUP FOR CHARGE 
SCATTERING 

In this section, we set up the numerical formalisms for 
charge scattering in /(i?) gravity, including field equa¬ 
tions, initial conditions, boundary conditions, discretiza¬ 
tion scheme, and tests of numerical codes. 


D. /(R) model 

For a viable dark energy f{R) model, /' has to be 
positive to avoid ghosts [54] , and /" has to be positive to 
avoid the Dolgov-Kawasaki instability |^. The model 
should also be able to generate a cosmological evolution 
compatible with the observations |56l HZ] and to pass 
the Solar System tests [39l . Equivalently, general 

relativity should be restored at high curvature scale, and 
the f{R) model mainly deviates from general relativity 
at low curvature scale comparable to the cosmological 
constant. In this paper, we take a typical dark energy 
f{R) model, the Hu-Sawicki model, as an example. This 
model reads |MI 


f{R) =R-Rc 




D2R^ + R^ 


(29) 


where n is a positive parameter, Di and D 2 are dimen¬ 
sionless parameters, Rq = 87rpo/3, and po is the average 
matter density of the current Universe. We consider one 
of the simplest versions of this model, i.e., n = 1, 


f{R) = R- 


DRaR 


R Rq 

where is a dimensionless parameter. In this model, 


(30) 


/' = 1 - ^^0 


{R + Ro?’ 


(31) 


A. Field equations 


In double-null coordinates (26), using 


Gl^G% = Stt 


^(tot3-l)i ^ ^(tot3,l)3 
one obtains the equation of motion for r, 


r{-r,tt + r,xx) - r'ft + = e ^1 - , 

where = dr/dt and other quantities are defined ana¬ 
logously. For simplicity, we define rj = and integrate 
the equation of motion for 77, instead m- The equation 
of motion for rj can be obtained by rewriting Eq. (331 as 


- V,tt + V,xx = 2e-2" ( I - STrrV - ^ ) . (34) 


Gg = provides the equation of motion for tr. 


~ cr -I- a^xx + 


' ,tt ' ,XX 



(35) 


In double-null coordinates, the equations of motion for 
(j) {21) and tp {22) become, respectively. 


V'{P>) = 


i?3 


Venf'^iR + Roy 




(32) 

As implied in Eq. (321, to make sure that the de Sitter 
curvature, for which w ((/>) = 0, has a positive value, the 
parameter D needs to be greater than 1. In this paper, 
we set D to 1.2 and set Rn to 10“^ or 10“®. Then, 


together with Eqs. (231 and (321, these values imply that 
the radius of the de Sitter horizon is about a/I/Rq ~ 10^. 
Moreover, in the configuration of the initial conditions 
described in Secs. |IIIB| and |V^ the radii of the outer 
apparent horizons of the formed black holes are about 
2.1 and 3.7, respectively. [See Figs. if) and &)■] The 
potentials in the Jordan and Einstein frames are plotted 
in Figs. a) and|^b), respectively. 


where 


“t“ 4^,xx H" ( 
r 


= e 


-2a- 


V'{cP) + 


(36) 


- 4’,tt + 'ip,xx + + r^x'tp,x) 

r 


(37) 




( 38 ) 


X 



























7 






Figure 3: Potentials for f{R) models. X=f' s-nd Xp=x/(1 + x)- U{x) and V(x) are the potentials in the Jordan and 

Einstein frames and can be obtained from Eqs. 0 and ( |23[ ), respectively.^ _ (a) and (b) are for the Hu-Sawicki model ( |30[ ), 

f{R) — R — DRoR/{R + Ro), while (c) and (d) are for the combined model ( |79[ ), f{R) ~ R — DRoR/{R + Ro) + aR^. D — 1.2, 
Rf) = 10“®, and a = 1. 


The {mm} and {mm} components of the Einstein equa¬ 
tions yield the constraint equations 




+ 2ct -I- 4nr 




(39) 


B. Initial conditions 

We set the initial data to be time symmetric: 

= cr,t = 4',t = ■0,* = 0 at t = 0. (43) 


-I- 2a,vr,v + 47rr 



(40) 


Via the definitions of m = {t — x)l2 and v = (t-|-a;)/2, the 
constraint equations can be expressed in (t, x) coordin¬ 
ates. Equations (40) — (39) and ( [40| ) -|- ([39| ) generate the 
constraint equations for the {te} and {tt} + {xx} com¬ 
ponents, respectively. 


Therefore, in this configuration, the constraint equa¬ 
tion (41) is satisfied identically. Note that, in this con¬ 


figuration, the values of r^t and a^t at t = 0 are the same 
as those in the Reissner-Nordstrom metric case. 

We set the initial value for ip as 


ip{x,t)\t=o = a ■ exp 


{x - Xof 


(44) 


r,tx + r,tcr,x + r,xCr^t + 47rM (t),t4>,x + 


X 


= 0, (41) 


r,tt + r,xx + ‘2'{r^tO-,t + r^xCr,x) 


i2 , 12 , t + 


+ Anr 4>, +4>,x + — 


= 0 . 


(42) 


In this paper, we give (p two sets of initial conditions: 

set 1 : (p{x,t)\t=o = (po, (45) 

set 2 : (p{x,t)\t=o = (po + , (46) 

where (po is de Sitter value, defined by V'{(po) = 0. The 



































initial value for a is defined to be the same as the corres¬ 
ponding one in the Reissner-Nordstrom case (27), 


it=o 


= e 


.- 2 f 7 |RN 

lt =0 


r+r_ nfc^r 


— - 1 
r_ 


l+T 


(47) 

where r is defined by Eq. (28) with t — Q. We obtain 
the initial value for r in charge scattering by combining 
Eqs. (ISSl) and (|42|), 


- + — 


2 r 


x2 , ,2 , 


- 27rr 4> , + 4>,x + — 


(48) 


2 r® 


1 


We set r^x = <y,x = ^ at the origin (x = 0,t = 0) as 
in the Reissner-Nordstrom metric case. The definition 
domain for the spatial coordinate x is [—Xb Xb\. Then 
r{x,t)\t=o can be obtained by integrating Eq. (48) via 


the fourth-order Runge-Kutta method from a; = 0 to 
X = ixb, respectively. The initial values of r, a, /', and 
^ are shown in Fig. |11| 

In this paper, we employ the finite difference method. 
The leapfrog integration scheme is implemented, which 
is a three-level scheme and requires initial data on two 
different time levels. With the initial data at t = 0, we 
compute the data at t = At with a second-order Taylor 
series expansion as done in Ref. [51]. Take the variable 
Ip as an example, 


V’|t=At = ip\t=o + ^.t|t=oAt -h -■!/'.tt|t=o(At)^. (49) 


The values of 4>\t=o and V',t|t=o are set up as discussed 
above, and the value of 'p, tt\t= o can be obtained from the 
equation of motion for ip (37). 


Up to this point, the initial conditions are hxed, with 
all the field equations being taken into account. The hrst- 
order time derivatives of r, u, p, and ipatt = 0 described 
by Eq. (43) ensure that the constraint equation (41) is 


satisfied. The equation for r^xx at t = 0 expressed by 
(48) implies that the constraint equation (42) is satisfied. 


Computations of r, a, p, and p at t = At via a second- 
order Taylor series expansion, as expressed by Eq. (49) 
for the case of p, satisfy all the equations of motion. 

Note that the region of a: < 0 is included in the initial 
conditions. This may not be physical. However, we are 
mainly interested in the interior dynamics of black holes, 
and then it is not important where the scalar field origin¬ 
ally comes from. This setup makes it convenient for us 
to compare the results of charge scattering to the known 
solutions of the Reissner-Nordstrom geometry. There¬ 
fore, we use this setup as a toy model. 


A 





^--.j 



I_L_i_J_! 


I_I 

Ax 


Figure 4: Numerical evolution scheme. 


C. Boundary conditions 

The values of r, a, p, and p at the boundaries of x = 
ixf, are obtained via extrapolations. In fact, since we 
are mainly concerned with the dynamics around x = 0, 
the boundary conditions will not affect the dynamics in 
this region, as long as Xb is large enough. 


D. Discretization scheme 


In this paper, we implement the leapfrog integration 
scheme, which is second-order accurate and nondissipat- 
ive. We let the temporal and spatial grid spacings be 
equal, At = Ax. 

The equations of motion for p {36 \ and p ( |T^ are 
coupled. Newton’s iteration method is employed to 
solve this problem [51]. With the illustration of Fig. 
the initial conditions provide the data at the levels of 
“down” and “here”. We need to obtain the data on the 
level of “up”. We take the values at the level of “here” to 
be the initial guess for the level of “up”. Then, taking 
p as an example, we update the values at the level of 
“up” using the following iteration: 

.xnew _ , _ G{pup) 

^up -Uup 

where G{pup) is the residual of the differential equation 
for the function pup, and J{pup) is the Jacobian defined 
by 


J (V’up) = 


dG{Pu 

dp- 


up 


We do the iterations for the coupled equations one by one, 
and run the iteration loops until the desired accuracies 
are achieved. 
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Figure 5: Tests of numerical code for charge scattering, (a) Numerical vs analytic results for a Reissner-Nordstrom black hole, 
m = 1, g = 0.7, and Ax = At = 10“^. The slice is for {x = 3Ax,t = t). This is a special case of charge scattering with 
contributions of scalar fields being set to zero. Numerical and analytic results match well at an early stage, while at a later 
stage gravity and electric held become stronger, the numerical evolutions have a time delay, compared to analytic solutions, 
(b) Numerical tests for the {tx} constraint equation (411 and the evolution of ijj on the slice {x = x,t = 0.65). They are both 
second-order convergent. 


E. Tests of numerical code 


To make sure that the numerical results are trust¬ 
worthy, one needs to test the numerical code. We com¬ 
pare the numerical results obtained by the code with the 
analytic ones for the dynamics in a Reissner-Nordstrom 
geometry, and examine the convergence of the constraint 
equations and dynamical equations in charge scattering. 


The dynamics in a Reissner-Nordstrom geometry is a 
special case for charge scattering, in which the contribu¬ 
tions from the scalar fields are set to zero. This special 
case has analytic solutions expressed by Eqs. (27) and 
(28). Therefore, we can test our code by comparing the 
numerical and analytic results in the Reissner-Nordstrom 
geometry. Set m = 1, q = 0.7, and Aa: = At = 

We plot the evolutions of r and tr along the slice {x = 
3xl0“^,t = t) in Fig. [^a). As shown in Fig. [^a), nu¬ 
merical and analytic results match well at an early stage; 
while at a late stage where gravity and electric field be¬ 
come strong, the numerical evolutions have a time delay 
compared to the analytic solutions. 


When the numerical results are obtained, we substi¬ 
tute the numerical results into the discretized equations 
of motion and the constraint equations, and find that 
they are well satisfied. See Figs. [I^ and [TS] for example. 
Moreover, the convergence of the constraint equations 
(41) and (42) is examined. We assume one constraint 
equation is nth-order convergent: residual=(!l(h"), where 
h is the grid size. Therefore, the convergence rate of the 
discretized constraint equations can be obtained from the 


ratio between residuals with two different step sizes, 


n = log 2 


0{h^‘ 


OKI) 


(50) 


Our numerical results show that both of the constraint 
equations are about second-order convergent. As a rep¬ 
resentative, we plot the results for the {tx} constraint 
equation (41) in Fig. [^b) for the slice {x = x,t = 0.65). 

Convergence tests via simulations with different grid 
sizes are also implemented [651 166] . If the numerical 
solution converges, the relation between the numerical 
solution and the real one can be expressed by 


^real = ^'’ 


0 (h"), 


where is the numerical solution. Then, for step sizes 
equal to /i/2 and /i/4, we have 


^real 


= F'i +0 



'^real 


= F-^ +0 



Defining ci = F^ — F 5 and C 2 = F 2 — F 4 , one obtains 
the convergence rate 


n = logs ■ (51) 

The convergence tests for r] = ^ a, (j), and ip are in¬ 

vestigated. They are all second-order convergent. As a 
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representative, the results for ip are plotted in Fig. [^b) 
for the slice {x = x,t = 0.65). The values of the 
parameters in charge scattering in this section are de¬ 
scribed at the beginning of Sec. lYll We use the spa¬ 
tial range of a: S [—10 10 ] and the grid spacings of 
h = Ax = At = 0 . 02 . 


IV. NEUTRAL SCALAR COLLAPSE 


In the numerical integration, once the value of r at the 
advanced level is obtained, the value of a at the current 


level will be computed using Eq. (57). 
We set the initial data as 


r,tt = = 4>,t = ' 0 ,t = 0 at t = 0 . 


(59) 


The initial values for %[= exp(-^2/3K0)] and ii){r) are 
defined as 


In this section, we consider neutral collapse in flat 
geometry in f{R) gravity and discuss the mass inflation 
which happens in the vicinity of the central singularity 
of the formed black hole. 


A. Numerical setup 

The numerical setup in neutral scalar collapse in /(i?) 
gravity is discussed in Ref. m- The dynamical equations 
for r, 77 , a, (p, and ip can be obtained by setting the terms 
related to the electric field in the corresponding equations 
in Sec. IIII Al to zero: 

r{-r,tt + r^xx) - r'^t + - Sirr'^V), (52) 


-V,tt+V,xx = ‘2e ^'^(l - STrr^V), (53) 


~ cr.tt + a XX + 


Xm - r. 



(54) 


x(r)|t=o = a • [I - tanh(r - rif] + xo, (60) 


7 /'(r)|t=o = & ■ tanh(r - r 2 )^, (61) 

with a = 0 . 2 , b = 0 . 1 , ri = r 2 = 4, and U'{xo) = 0 . 
The parameters for the Hu-Sawicki model (30) are set as 
D = 1.2 and Rq = 10"®. 


The local Misner-Sharp mass m is defined as m 


2m 




(62) 


(See Ref. [55] for details on various properties of the 
Misner-Sharp mass/energy in spherical symmetry.) Then 
on the initial slice {x = x,t = 0 ), the equations for r, m, 
and g are [201 ED 


r,a: = ( 1 - e® 


m r = 47 rr^ 



(63) 


(64) 


~ 4‘,tt + (p,xx H- { — X^t(p,t + X^x(p,x) 


= e 


-2(T 


v'{(p) + -pKT 

v 6 


- i’,u + '(p,xx + -{-r,t'fp,t + r,xtp,x) 
r 


= \l ■^K[-(p,ti’,t+ ^,xtp,x), 


(55) 


(56) 


where = e^'^{ip‘^^ — ip'^x)/x- 

In the equation of motion for a ( [54] ), the term {r^tt — 
x,xx)lx can create big errors near the center x = r = 0. 
To avoid such a problem, we use the constraint equation 


(39) alternatively m- Defining a new variable g 

g =-2 ct - ln(-r,„), (57) 


one can rewrite Eq. (39) as the equation of motion for g, 

ip‘^ \ 

' (58) 


= dTT • — • [(/)_„+ ^ 


5,1 



= 47rr (px H-^ 


(65) 


Set r = m = g = 0at the origin {x = 0,t = 0). Then the 
values of r, m, and g on the initial slice (x = x,t = 0) can 
be obtained by integrating Eqs. ([63|)-(65) from the center 


a: = 0 to the outer boundary x = xt via the fourth-order 
Runge-Kutta method. The values of r, tr, </>, and ip at 
t = At can be obtained with a second-order Taylor series 
expansion, as discussed in Sec. [nra The value of g at 
t = At can be obtained using Eq. (57). 


The range for the spatial coordinate is defined to be 
a: G [0 20]. At the inner boundary a: = 0, r is always set 
to zero. The terms 2{—r^t4‘,t + x x<P,x)/x in Eq. (55) and 
2{—r tip t -\- r^x'>P,x)/x in Eq. (56) need to be regular at 
X = r = 0. Since r is always set to zero at the center, so 
is r^f Then we enforce (p and ip to satisfy (p^x = '<P,x = 0 
at X = 0. The value of g at x = 0 is obtained via ex¬ 
trapolation. We set up the outer boundary conditions at 
X = 20 via extrapolation. The temporal and the spatial 
grid spacings are At = Ax = 0.005. 

The numerical code is second-order convergent and is 
the one developed in Ref. [50] . 
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Figure 6: Evolutions in neutral collapse for the Hu-Sawicki model (301, f{R) = R — DRqR/(R + Rq). In (a)-(d), the time 
interval between two consecutive slices is lOAt = 0.05. (e) and (f) are for the apparent horizon and the singularity curve of the 
formed black hole. 


B. Black hole formation 


The evolutions of r, cr, </>, and ip are plotted in 
Figs. [la)-[§d), respectively. On the apparent horizon of 


a black hole, the expansion of the outgoing null geodesics 
orthogonal to the apparent horizon is zero |69j . Then in 
double-null coordinates, on the apparent horizon, there 


















12 




b 

;h 

,o 


o 

W 




10 ^° 

10 ® 

10 ® 

lO'* 

10 ^ 

10 ® 



- - o-.ti 

V> 

1 ^,XX \ 

---\r,u/r\ 


. 1 - 


— 


--- 


llTTb^/'l 



__ 

-87re-2‘^F| 


- 10 ® residual 




0 0.01 0.02 0.03 0.04 


r 

(c) 




Figure 7: (color online). Dynamics on the slice (x = l,t = t) in neutral collapse for the Hu-Sawicki model (30l. (a)-(c); 
dynamical equations for r, t], and cr. (b) Near the central region, due to accumulation, the scalar field (j) is strong, a is positive. 
As a result, in the equation of motion for 77 (531, the terms 26“^°^ and 16TTe~^°^r^V are negligible. The equation is reduced 
to T]^u ~ ri,xx- (d) In(mEF) = + b, a = —1.6438 ± 0.0008, b = —0.998 ± 0.003. In(mjF) = alnr + b, a = —2.385 ± 0.002, 

b = -3.443 ±0.009. 


is [TD] 

± = 1 -^ 

Using this property, we locate the apparent horizon and 
plot it in Figs. [^e) and if)- As shown in Fig. ia), 
the central singularity is also approached in the collapse. 
Therefore, a black hole is formed. 


C. Asymptotic dynamics in the vicinity of the 
central singularity of the formed black hole 

We focus on the dynamics in the vicinity of the central 
singularity of the formed black hole. We will discuss that, 
in the vicinity of the singularity, due to the backreaction 
of the scalar fields on the geometry, the Misner-Sharp 
mass diverges. In other words, in addition to the in¬ 
ner horizons of Reissner-Nordstrom and Kerr black holes. 


mass inflation also happens in the vicinity of the central 
singularity of a Schwarzschild black hole. 

In the vicinity of the central singularity, the field equa¬ 
tion can be reduced to the following forms [ID] : 


rr^tt R 


(67) 

CTtt R 

r 

(68) 

4>,tt « 

. 2 , 

r 

(69) 

V'.tt « 

2 12 

(70) 
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(c) 
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Figure 8 : (color online). Dynamics on the slice {x = 2,t = t) in neutral collapse for the Hu-Sawicki model, (a)-(c): dynamical 
equat ions for r, 77 , and a. (b) At large-® regions, the scalar field (p is weak, and a is negative. As a result, in the equation of motion 
for 77 ( 53 1 , the term 26“^“^ is important. The equation is reduced to 77,41 ~ — 26“^°^. (d) In 771 = alnr-|- 6 , a = —0.2673 ± 0.0008, 
6 = 0.217 ±0.003. 


The asymptotic solutions to Eqs. (67)-(69) are 
r « 


(71) 


(T « S In^ ± CTo « [/3(1 — /?) — dTrC^] In^ ± ctq, (72) 


Cln^. 


(73) 


The variable ^ is defined as ^=to ~ where to is the 
coordinate time on the singularity curve. 

As implied in Eq. (70), ip is suppressed (magnified) 


by (p when cp^t is negative (positive). Due to the com¬ 
plex competition between gravity and dark energy, an 
approximate analytic expression for ip is not obtained. 
Substituting Eq. 0 into Eq. ( [^ yields 


Then we have 




(74) 


Using Eqs. 0 and ( [ 7 ^ , there 
X « 


IS 


(75) 


So as shown in Fig. |^c), when C is positive, x ^-Iso 
approaches zero as ^ and r approach zero. Then the 
transformation between the Jordan and Einstein frames, 
breaks down. Actually this is not a seri¬ 
ous problem, because numerical simulation stops anyway 
when the central singularity is approached. Moreover, in 
this paper, we discuss the dynamics in the vicinity of the 
central singularity rather than on the central singularity. 


D. Mass inflation 


In the vicinity of the singularity, in the equation of mo¬ 
tion for a ( 68 ), because of the contribution from cp, a{x, t) 
is greater than the corresponding value in the Schwarz- 
schild black hole case. This makes the mass function 
divergent near the singularity as will be discussed below. 
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Near the singularity, using Eqs. 
mass function can be written as 


.= [I+ 


-(l-if2)A3g2<To ^ 




, @-([7|), the 


(76) 


I _ ^2^^3+167rC^g2CTo 
8 


III 

& 

r,x 


4^,x 





JF 

9,t 


where K = \r^x/'i',t\- In the Schwarzschild black hole 
case, (7 = 0. The mass function is always constant and 
is equal to the black hole mass. In neutral collapse, the 
parameter /? does not change much and remains about 
1/2. However, the parameter C is not zero. Then the 
metric quantity a is modified. [See Eq. (72).] As a res¬ 
ult, the delicate balance between r and + r\) is 


broken. Consequently, as implied in Eq. (76), near the 
singularity, the mass function diverges: mass inflation 
occurs. 

During the collapse, before the black hole is formed, 
the energy of the scalar fields accumulates in the cent¬ 
ral region. As a result, the scalar fields near a; = 0 are 
stronger than those at large-a: regions. Next we discuss 
three consequences. As a support, we examine the dy¬ 
namics in the vicinity of the singularity via mesh refine¬ 
ment that was implemented in Refs. uniin], and plot 
two sample sets of results on the slices {x = l,t = t) and 
{x = 2,t = t) in Figs. and 1^ respectively. 

(i) Values of a. Due to the backreaction of the scalar 
fields on the geometry, cr in small-a: regions is 
greater than in large-x regions. In fact, a is pos¬ 
itive in small-x regions, while negative in large-x 
regions. [See Fig. [^b).] Our numerical results 
of the parameter C in (j>KiC\n^ are C « 0.18 at 
X = 1 and C « 0.07 at x = 2. Then we have 

47r(7^ « 0.41 > 1/4 and a > 0 at x = 1, while 

« 0.06 < 1/4 and cr < 0 at x = 2. [See 


Eq. (72). 


(ii) Equation of motion for ij (53). For positive a, the 
terms 2e“^'^ and 167re“^'^rH4 in Eq. ( [5^ are neg¬ 
ligible, compared to the other two. Then Eq. (53) 
is reduced to r]^tt ~ ri,xx- [See Fig. [^b).] How¬ 
ever, for nega tive a, the term 26“^*^ is important, 
and Eq. ( [^ is reduced to i] u ~ —26“^*^. [See 
Fig-ib).] 

(iii) Growth of the mass function. As implied in 
Eq. ([7§), the mass function grows faster in the 
strong scalar field case than in the weak one. [See 
Figs. [3d) andjfd).] 

Since f{R) gravity is defined in the Jordan frame, it is 
interesting to examine the mass function in the Jordan 
frame. A generalized Misner-Sharp energy in f{R) grav¬ 
ity in the Jordan frame was defined in Ref. [40) . However, 
due to the complexity of some integrals, an explicit quasi¬ 
local form is usually not available with the exceptions of a 


Friedmann-Robertson-Walker universe and static spher¬ 
ically symmetric solutions with constant scalar curvature. 
In this paper, for simplicity, we remain to use the conven¬ 
tional format of the Misner-Sharp function. Considering 
the transformation that we used, there 

are e^'^jjp = X’S^'^Ief and rjp = Considering 

that near the central singularity uniiii] 


(77) 


one can obtain that iXjpwAlEF- Then the mass function 
in the Jordan frame can be written as 

wjF = ^ [1 + - ?'^^)|jf] 

" 2 “ 



(1 - 




1^1 _ (1 - 


fKCg2cTo 


(78) 

In the case of (7 > 0, due to the factor ^ mjp is 

greater than toef- The tojf along the slice (x = l,t = t) 
is plotted in Fig. Hd). 

For stationary black holes (e.g., Schwarzschild and 
Reissner-Nordstrom), the Misner-Sharp mass function is 
always equal to the black hole mass. In spherical sym¬ 
metry, at spatial infinity, the mass function describes 
the total energy/mass of an asymptotically flat space- 
time [68) . In gravitational collapse case, it means the 
total mass of the collapsing system. 

In the vicinities of the central singularity of a Schwar¬ 
zschild black hole and the inner horizon of a Reissner- 
Nordstrom or Kerr black hole, the dynamics and some 
quantities are local. The mass function is just a para¬ 
meter which varies at each point, not giving global in¬ 
formation on the black hole mass. 


V. NEUTRAL SCALAR SCATTERING 

In this section, we consider neutral scattering, in which 
a neutral scalar field collapses in a Schwarzschild geo¬ 
metry in f{R) gravity. The numerical formalism is a 
simpler version of the one in charge scattering that has 
been constructed in Sec. jlll) and it can be obtained by 
removing the electric terms in the field equations presen¬ 
ted in Sec. [Ill A) and replacing the Reissner-Nordstrom 
geometry with a Schwarzschild one. 


A. A dark energy f{R) singularity problem 

In neutral scattering, for usual initial conditions of the 
scalar degree of freedom /', f' asymptotes to zero as 
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the central singularity is approached, which is similar to 
what happens in the neutral scalar collapse discussed in 
Sec. |IV| Details are skipped here. On the other hand, 
when the initial velocity or acceleration of /' is large 
enough, /' can become 1 before the central singularity 
is approached. As implied in Eq. ( [^ , for dark energy 
f{R) models, this means that the Ricci scalar becomes 
infinite, and the simulation breaks down. In this section, 
we will focus on this singular circumstance. 

The parameters are set as follows: 

(i) Schwarzschild geometry: m = 1. 

(ii) Physical scalar field: = a ■ 

exp [—(a: — Xo)‘^/b\ , a = 0.08, 6=1, and xq = 4. 

(hi) f{R) model: f{R) = R-DRoR/{R +Rq), D = 1.2, 
and Rq = 10“®. 

(iv) Scalar degree of freedom: (j){x,t)\t=o = a + b- 
exp [—(cc — a:o)^/c], a = —0.02, 6 = —0.1, c = 1, 
and Xq = —2. 

(v) Grid. Spatial range: x € [—10 10]. Grid spacings: 
Ax = At = 0.005. 

The results in this circumstance are plotted in Fig. 
We plot the dynamics of </> on the slice (x = —2, t = t) in 
Fig.^d), from which one can see that (j) is mainly accel¬ 
erated by the spatial derivative and the geometrical 
term —2r^t<t’,t/x. Eventually, (j) and f' approach 0 and 1, 
respectively. Then the Ricci scalar becomes singularity, 
and the simulation stops, as shown in Fig.j^a). 


B. Avoidance of the singularity problem 


In fact, it has been argued that such a singularity 
problem can also be caused in cosmology and compact 
stars m- This problem can be avoided by adding an 
R^ term to the dark energy f{R) model [1^49j . In this 
paper, we add the R^ term to the Hu-Sawicki model, 

f{R) = R-^^+aR\ (79) 

JrC -h Kq 


In this combined model, at high curvature scale, f' « 
2aR. So a singular R is pushed to regions where /' 
and (j) are also singular. Then the singularity problem 
is avoided. The parameters take the same values as in 
the last subsection. In addition, the new parameter a in 
the /(i?) model (791 is set to 1. 

The numerical results for neutral scattering for this 
modified model are plotted in Fig. As shown in 

Fig. [^b), for this model, f can cross 1 without dif¬ 
ficulty. What are the limits that f and R can reach? 
As shown in Fig. [l0j[d), as the central singularity is ap¬ 
proached, there is 


4>,tt « 



(80) 


Since r^t is negative, the above equation describes a pos¬ 
itive feedback system of (p^tt and \4>,t\ produces more 
of and in turn produces more of \4>^t\- Since 

4>^t is positive, p can be rapidly accelerated to positive 
infinity. Then /' and R will go to infinity as the central 
singularity is approached. This feature surely deserves 
some attention as discussed below. 

The Starobinsky model, f{R) = R + aR?, was ob¬ 
tained by taking into account quantum-gravitational ef¬ 
fects. It could cause inflation in the early Universe and 
is conventionally believed to be singularity-free [H]. The 
combined model is reduced to the Starobinsky model at 
high curvature scale. It is found that, in neutral scatter¬ 
ing for the combined model, a new black hole, including 
a new central singularity, can be formed. Near the cent¬ 
ral singularity, gravity dominates other terms, including 
the potential related to the R? term, such that the Ricci 
scalar R can be pushed to infinity by gravity. We also 
simulate scalar collapse for the Starobinsky model in flat 
geometry. Similar results are obtained m- Therefore, 
the classical singularity problem, which is present in gen¬ 
eral relativity, remains in collapse for these models. Fur¬ 
ther details are skipped. 


VI. RESULTS FOR CHARGE SCATTERING 


In this section, we explore charge scattering: scalar 
collapse in a Reissner-Nordstrom geometry. We study 
the evolutions of the metric components and scalar fields 
and obtain approximate analytic solutions. We closely 
compare the dynamics in Schwarzschild black holes, 
Reissner-Nordstrom black holes, neutral scalar collapse, 
and charge scattering. 

In this section, the parameters are set as follows: 

(i) Reissner-Nordstrom geometry: m = 1, and q = 0.7. 

(ii) Physical scalar field: tp{x,t)\t=Q = a ■ 

exp [—(x — XQ)'^/b] , a = 0.08, 6=1, and xq = 4. 

(hi) fiR) model: f{R) = R- DRqR/{R + Rq), D = 1.2, 
and Rq — 10“®. 


(iv) Scalar degree of freedom: (j){x,t)\t=o = 4 >q, with 
v’{Pq) = 0 . 

(v) Grid. Spatial range: x S [—10 10]. Grid spacings: 

VI A| and VI B[ and 


Ax = At = 0.005 for Sec s. 

Ax = At = 0.0025 for Secs. IVI Gl and IVI PI 


A. Evolutions 


1. Outline 


In this subsection, we describe the evolutions of r, tr, 


/', and Ip that are plotted in Fig. 11 Examining the 
equations of motion (33)-(37) and the numerical results 
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Figure 9: (color online). Results for neutral scalar collapse in a Schwarzschild geometry for the Hu-Sawicki model, (a) and (b): 
evolutions of r and /'. The time interval between two consecutive slices is 30At = 0.15. (c) evolutions of r and /' on the slice 
(x = —2,t = t). (d) dynamical equation for <j) on the slice {x = —2,t = t). 4> is mainly accelerated by the spatial derivative 
and the geometrical term —2r^t(p,t/r. Eventually, (p and f' approach 0 and 1, respectively. Then the Ricci scalar becomes 
singularity, and the simulation stops. 


plotted in Figs. [T^ and [15] one can see that in the charge 
scattering dynamical system, there are three types of 
quantities as follows: 


(i) Metric components: r and a. 
gravity. 


They contribute as 


(ii) Scalar fields: (j) and ip. 
gravitating fields. 


They contribute as self- 


(iii) Electric field and V{(j)). As implied in Eq. (331, 
they are repulsive forces. However, the numerical 
results show that in charge scattering, compared to 
the contributions from other quantities, the contri¬ 
bution from V{(j)) is negligible. 


Furthermore, these quantities can be separated into two 
sides: the gravitating side (r, a, (p, and ip) and the repuls¬ 
ive side [electric held and V{(p)]. The dynamics in charge 


scattering consists mainly of how these variables interact 
and how the gravitating and repulsive sides compete. 

According to the strength of the scalar held, charge 
scattering can be classihed into hve types as follows: 


(i) Type I: spacelike scattering. When the scalar held 
is very strong, the inner horizon can contract to zero 
volume rapidly, and the central singularity becomes 


spacelike. Sample slice: {x = 1.5, f = t) in Fig. 11 
See Sec. IVIBI 


(ii) Type II: null scattering. When the scalar held is in¬ 
termediate, the inner horizon can contract to a place 
close to the center or reach the center. For each vari¬ 
able, the spatial and temporal derivatives are almost 
equal. In the case of the center being reached, the 
central singularity becomes null. This type has two 
stages: early/slow and late/fast. In the early stage, 
the inner horizon contracts slowly, and the scalar 
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Figure 10: (color online). Results for neutral scalar collapse in a Schwarzschild geometry in a combined f{R) model (791, (a) 
and (b): evolutions of r and f'. The time interval between two consecutive slices is 30Af = 0.15. (c) evolutions of r and f' on 
the slice {x = —2, t = t). (d) dynamical equation for (j) on the slice (x = —2, t = t). 


field also varies slowly. In the late stage, the in¬ 
ner horizon contracts quickly, and the dynamics is 
similar to that in the spacelike case. Sample slice: 


(x = 0.5, t = t) in Fig. 11 See Secs. VIC and VID 


In this paper, we will discuss Types I, II, and IV. 


2. Causes of mass inflation and evolutions 


(iii) Type III: critical scattering. This case is on the 
edge between the above two cases. When the central 
singularity is reached, it becomes null. Sample slice: 
(x = 1.4225, t = t) in Fig. 


11 


Due to the similarity 
to that in general relativity discussed in Ref. [41], 
details on this type of scattering in /(i?) gravity are 
skipped in this paper. 


(iv) Type IV: weak scattering. When the scalar field 
is very weak, the inner horizon does not contract 
much. Sample case: Fig. See Sec. |VII| 

(v) Type V: tiny scattering. When the scalar field is 
very tiny, the influence of the scalar field on the 
geometry is negligible. Sample slice: (x = —3, t = t) 
in Fig. 


The local Misner-Sharp mass for a charge black hole 
hole is 


5 " 


,/i ' ,1/ 


= = + (81) 


In a Reissner-Nordstrom geometry, in Kruskal-like co¬ 
ordinates expressed by Eq. (27), near the inner horizon, 
although a asymptotes to positive infinity, — r\) is 


much less than Consequently, — rt^) ap¬ 

proaches zero. Then as implied in Eq. (pR, the mass 
function m takes a hnite value and is equal to the black 
hole mass. For more details, see Ref. [4l] . 

In charge scattering, the equations of motion for ri=r^ 
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X 

(a) 



X 

(c) 



X 

(e) 



X 

(b) 



X 

(d) 



X 

(f) 


Figure 11: Evolutions in charge scattering, (a)-(d): evolutions of r, a, /', and ip. The time interval between two consecutive 
slices is 30At = 0.15. (e) and (f) are for the apparent horizon and the singularity curve of the black hole. When the scalar field 
is strong enough (around x = 2), the inner horizon can be pushed to the center, and the central singularity becomes spacelike. 
When the scalar field is weak enough (e.g., —4 < x < —1), the inner horizon does not change much. At the intermediate state 
(e.g., 0 < a; < 1.5), the inner horizon contracts to zero, and the central singularity becomes null. The results for the inner 
horizon, especially for x > 2, are not that accurate. We are aware that the inner horizon is actually at infinity, while r still can 
be very close to r_ when x and t take moderate values. 
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and a are 


- ??,tt + V,xx = 2e ( 1 - Sirr'^V - ^ y , 


, , T,it - r 

^,tt + cr XX H - 

r 

/,2 _ J ,2 


(82) 


■ 47r 
.-2 <t9 


X 


- 2e-^'"y 


2 _ ^ 2 ^ ^ i’,x 
2 


(83) 


e = 0. 

J.4 


At the beginning, \(j)^t\ and |'(/',t| naay be less than \(j)^x\ 
and respectively. However, as r decreases toward 

the central singularity, gravity becomes stronger. Then 
and become greater than \4>^x\ and \ip,x\, re¬ 
spectively. [See Fig. 15 ] As a result, in this case, the re¬ 
pulsive “force”, 47r[<()^^ (j)% + (V'^t ~ V'^a;)/x] + 
is greater than the corresponding one, in a 

Reissner-Nordstrom geometry. This makes a accelerate 
faster than in the Reissner-Nordstrom geometry. Con¬ 
sequently, the repulsive force from — 1) for 

r 7 =r^ is much weaker than the corresponding value in the 
Reissner-Nordstrom geometry. As a result, near the inner 
horizon, jr^tj is much greater than the corresponding one 
in the Reissner-Nordstrom case. Then {r\ — r\) moves 
from extremely tiny values in the Reissner-Nordstrom 
metric case to moderate values, and r crosses the inner 
horizon r = r- for the given Reissner-Nordstrom geo¬ 
metry. With Eq. (81), the mass parameter grows dra¬ 
matically: mass inflation takes place. In other words, 
regarding the causes of mass inflation in charge scatter¬ 
ing, the scalar fields’ backreaction on r is more important 
than that on a. The evolutions of r, a, /', and '0 are 
plotted in Figs. [II|[a)-|TI|[d). 

Because of gravity from the black hole and the contri¬ 
bution from the physical scalar field, for some configur¬ 
ations, /' can decrease to zero as the central singularity 
is approached. [See Eq. (36) and Eigs. 11 13 and 15] 


Nothing is wrong with this circumstance. However, for 
certain configurations, f' can be pushed to 1, as shown 
in Eig. Correspondingly, the Ricci scalar R becomes 
singular. We will discuss the former case in this section 
and the latter case in Sec. Enn 

The evolution of tp is plotted in Fig. [TI|(d). In the 
configurations that we consider, as the central singularity 
is approached, because of the strong suppression from the 
dark energy scalar 0, ip asymptotes to constant values. 


3. Locations of horizons 

We locate the outer and inner horizons using the fol¬ 
lowing equation, 

+ r^) = I- —+= Q. (84) 


The results are plotted in Figs. 11 [e) and [T^f). Due 
to the absorptions of the energies of (p and ip, the outer 
horizon increases from the original value of 1.7 to 3.7. 
Note that the results for the inner horizon, especially at 
regions where x > 2, are not that accurate. We are aware 
that the inner horizon is actually at infinity, while r still 
can be very close to r_ even when x and t take moderate 
values. 


B. Spacelike scattering 

In spacelike scattering, the scalar field is so strong, 
such that the inner horizon can contract to zero volume 
rapidly, and the central singularity converts from timelike 
into spacelike. Taking the slice {x = 1.52, t = t) as an 
example, we plot the evolutions of r, a, f, ip, and m on 
this slice in Fig. The mass function remains equal to 
the mass of the original Reissner-Nordstrom black hole, 
mo = 1, until r is very close to r_. By then mass inflation 
takes place. We examine the dynamics near the central 
singularity via mesh refinement and plot the terms in the 
dynamical equations in Fig. [TSl 

The strongness of the scalar field causes several con¬ 
sequences as below. 

(i) Motion of r. The quantity r does not decelerate 
much when it crosses the inner horizon of the given 
Reissner-Nordstrom black hole, and it can approach 


the center. [See Fig. 12 


and the term e 


(ii) Nature of the central singularity. The central sin¬ 
gularity converts from timelike into spacelike. 

(hi) The dynamics in spacelike scattering is similar to 
that in strong, neutral scalar collapse. The quant¬ 
ity a takes large positive values, such that in the vi¬ 
cinity of the central singularity, compared to o ther 
’’ jr"^ in the equation for r (33) 
^q^ jr^ in the equation for cr (35) 
are negligible. As a result, in the vicinity of the 
central singularity, the dynamics is similar to that 
in strong, neutral scalar collapse as expressed by 
Eqs. (67)-(70). [See Figs. and 13] Then the 
quantities r, a, (p, and m take similar forms as those 
in neutral collapse. [See Fig, [m]] 

In both strong, neutral scalar collapse and spacelike 
scattering, the equation of motion for g is reduced 
to g^tt ~ 'n,xx- [See Figs. Bb) and [^b).] 

Since (p^t is negative, the term yj2j^K(p^t'4’,t in 
Eq. (37) functions as a friction force for ip. Con¬ 
sequently, compared to (p, ip grows slowly. As the 
singularity is approached, ip even approaches con¬ 
stant values. [See Figs, [l^c) and[l4j[d).] As shown 
in Fig. 13 e), near the singularity, two major sets of 
terms can be expressed as 

r,t'^P,t^r^xi’,x, (p,ti!,t~(p,xi’,x- 
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(a) (b) 

Figure 12: Evolutions along the slice {x = 1.52, t = t) in one spacelike scattering process, (a) evolutions of r and a. (b) 
evolutions of /', 'ip, and m. In this case, the scalar fields are so strong, such that r can rapidly cross the line of r = r_ and 
approach zero. Other variables {a, /', ip, and m) also evolve rapidly after r has crossed the line of r = r_. 


Alternatively, 


tA 

'^,x 


r,x 



(85) 


(iv) Growth of mass function. In spacelike scattering, 
the equation of motion for (p can be simplified as 

a^tt ^ ( 86 ) 

Consequently, with the results obtained in 
Sec. |IV C| a has the following asymptotic solution: 

a ~ B In^ + (To « —47rC^ In^ + (Tq, (87) 


with (/) « Clnf. Then similar to the mass function 
(76) in neutral collapse, with Eq. (81), the mass 
function in spacelike scattering can be written as 


r 

m = - 
2 




(1 - A:2)^3g2ao 




^(l-if2)^4(-B+l)g2.o 
1(1 _ ^2^^4(47rC= + l)g2<To 


„4B-1 


„-167rC"'-l 


( 88 ) 


Numerical results show that, for the sample slice 
(x = 1.52, t = t), near the central singularity, the 
slope of the singularity curve K is about 0.04. As 
shown in Fig |14K d), we linearly fit the numerical 
results of m via 

Inmwalnr + b, (89) 

obtaining 

a = -16.937 ± 0.003, 6 =-29.04 ± 0.01. 


Fitting numerical results for a according to Eq. (87) 
and combining Eqs. (88) and (89), we obtain 


flanalytic = 4i? - 1 = 17.0 ± 0.8, 


b 


analytic 


= In 


i(l-X2)yl4(-B + l)g2.o 


-28 ±4. 


Similarly, fitting numerical results for (p according 
to 0 « Cln^, we obtain 

Oanalytic = - 1 = 18.06 ± 0.01, 


b 


analytic 


= ln 


1(1 _ 1^2^^4(477C2 + l)g2<T0 
8 


-28 ±4. 


One can see that the above three sets of results 
match well. 


C. The late/fast stage of null scattering 

When the scalar fields are less strong, the inner horizon 
may still contract to zero. However, in this case, the 
central singularity becomes null rather than spacelike. 
The equations of motion remain null: they have similar 
forms as free wave equations, e.g., (p^tt ~ <P,xx- 

In the Reissner-Nordstrom black hole case, near the 
center, the repulsive (electric) force dominates gravity, 
and the central singularity is timelike. In spacelike scat¬ 
tering as discussed in the last subsection, gravity from 
the scalar field and the background geometry dominates 
the repulsive force. As a result, the central singularity 
is spacelike. At the late stage of null scattering that 
will be studied in this subsection, the scalar field is less 
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(e) 



(b) 



(d) 


Figure 13: (color online). Dynamics along the slice (x = 1.52, t = t) near the spacelike, central singularity. Near the singularity, 
we can approximately rewrite the original equations of motion for r, r/, a, (j), and ip as follows, (a) rr ,tt ~ -r'^f (b) ri,u ~ V%- 
(c) (T,tt « (d) (p.tt ~ -2r^t4>,tlr. (e) r^tip,t~r,xip,x and (j>,tip,t ~ (p.xip.x- Then we have ip,t/ip,x~r,x/r,t ~ 


strong, and the central singularity is null. Because of this, 
one may say that null scattering is a critical case of the 
competition between repulsive and gravitational forces, 
in which case the two types of forces have a balance. 

Take the slice {x = 0.5, t = t) as a sample slice, we 


plot the terms in the held equations for r, cr, <j), and ip in 
Fig. [l^ an d the evolutions of r, cr, cp, ip, m, and |1 — K'^\ 
in Fig. |16| We investigate the dynamics in the vicinity of 
the central singularity via mesh rehnement and plot the 
results in Fig. 
























































22 









(a) 




Figure 14: Solutions along the slice (x = 1.52, t = t) near the spacelike, central singularity, (a) lnr~aln^+fe, a = 0.5032±0.0001, 
b = -0.3033 ± 0.0008. a^aln^ + b, a = -4.0 ± 0.2, b = -10 ± 2. (b) <?!>«aln^±6, a = 0.5827 ± 0.0002, b = -1.294 ± 0.002. 
lnm«alnr ± 6, a = —16.937 ± 0.003, b = —29.04 ± 0.01. (c) and (d): ip{x,t) asymptotes to constant values as the central 
singularity is approached. 


The null scattering has two stages: early/slow and 
late/fast. As shown in Figs. 15 and 16 at the begin¬ 


ning of charge scattering, because of the repulsive force 
from the electric field, r, a, (/), and tp evolve slowly. As 
a result, the mass function m also grows slowly. We call 
this stage the early/slow stage. Later on, as the center 
is approached, gravity becomes very strong. Then these 
quantities evolve faster. We call this stage the late/fast 
stage. 

As shown in Fig. 


natural to guess that the quantities r, a, ip, and m may 
have expressions similar to those in spacelike scattering. 
In fact, this guess is verified by the numerical results 
plotted in Fig. [T^ Then we have 


of motion for 
as 


15 when r is very small, the equations 


'■>,u 


r « 


(j ~ B In^ ± (Tg « —d-TrC^ In^ ± ctq, 


r (33l, cr 

(351, and ()> (36) can be rewritten 

(/>« ClnC, 

- rr^tt 

« -rr^xx^ 

2^2 

(90) 

m ^ 



XX « 

2 

2 

(91) 



~ 4^,XX ^ 

r 

~ '^,x4^,X’ 

r 

(92) 


^(1- 


(93) 


(94) 


(95) 


c2B-i 


„4B-1 


„-167rC^-l 


(96) 


Since the above three equations have some similarities 
to the corresponding ones in spacelike scattering, it is 


As shown in Fig. 16’f), at the late stage, (1 — K^) is 
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Figure 15: (color online). Dynamics along the slice {x = 0.5, t = t) in null scattering, (a)-(f): dynamical equations for r, t], a, 
cj>, and ij). 


around 10“"*^. Due to the similarity to spacelike scatter¬ 
ing that we discussed in the last subsection, a compar¬ 
ison between numerical and analytical results for the fast 
stage of null scattering is skipped. 


D. The early/slow stage of null scattering 


As shown in Fig. 15 


at the early/slow stage of null 
scattering, the equations of motion for r (331, a (351, (j) 
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(c) (d) 





Figure 16: Evolutions along the slice (x = 0.5,t = t) in null scattering, (a)-(f): evolutions for r, a, f, ip, m, and |1 — K‘^\. At 
the early stage of mass inflation, 1.2 <t<2,r varies slowly; while at the late stage, t > 2, r varies rapidly toward zero. 


(36), and ip [37) are reduced as follows: 


^,tt ~ ^,xxi f,xx 1 4*,x'i (9^) 




(97) 


(l>,xx, r^t4’,t^r^x4>,x; (99) 
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o numerical results i 
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(a) 






Figure 17: Solutions near the null central singularity along the slice {x = 0.5, t = t) in null scattering, (a) InrRialn^ + h, 
a = 0.4421 ± 0.0008, b = -0.887 ± 0.002. (j«aln^ + h, a = -5.15 ± 0.01, b = 27.26 ± 0.02. <?i«aln^ + 6, a = 0.583 ± 0.002, 
b = 0.465 ± 0.005. (b) iIj approaches a constant value 0.4. Inmssalnr + b, a = —25.2 ± 0.1, b = 17.9 ± 0.3. 


V',it ~ ~ 4’,xll’,x, ( 100 ) 

The above equations are like free scalar wave equations 
in flat spacetime. The derivatives of one variable (r, a, </>, 
and ijj) are independent from the derivatives of another. 
Arbitrary functions of {t + x) or (t — x) can satisfy the 
above equations, and in principle, the initial conditions 
right after the collision between the scalar fields and the 
inner horizon will decide which function each variable 
can take. On the other hand, we find that, as shown in 
Fig. 18 the constraint equations (41) and (42) provide 


some useful information on the connections between some 
variables at the early stage of charge scattering: 


~ 47rr_^. 


( 101 ) 


As shown in Fig. 16 the quantities r, cr, /', ip, |1 — 


and m change dramatically at the beginning of charge 
scattering where r«r_. Note that, near the central sin¬ 
gularity, r, a, and /' have approximate analytic expres¬ 
sions in terms of ^ where to is the time coordinate 

of the singularity curve. So it is natural to guess that, 
at the early stage of charge scattering, the above quant¬ 
ities may also have approximate analytic expressions of 
Q = t — ts, where tg is a certain time value related to the 
early stage of charge scattering. We plot the evolutions 
of these quantities at the early stage of charge scattering 
in Fig. |18[ from which one can see that r, r^t, and a may 
have the following approximate analytic expressions: 


r,t-aC , 

(T«61nC -I- (To- 


( 102 ) 

(103) 

(104) 


In fact, a logarithmic expression for a is supported 
by its behavior near the inner horizon in the Reissner- 
Nordstrom black hole case. From Eqs. (28) and (47), one 


obtains that, as r approaches the inner horizon r = r_, in 
the case of t^x, a can be approximated by a logarithmic 
function of t. As shown in Fig.[T^c), describing the terms 
in the equation of motion for a in charge scattering, at 
the very early stage (t « 1) of the collision between the 
scalar fields and the inner horizon, compared to those 
from other terms, the contributions from the terms re¬ 
lated to (p and Ip are tiny. Therefore, at this stage, the 
evolution of a should not be much different from the cor¬ 
responding one in the Reissner-Nordstrom geometry. 

As shown in Figs.[I^c) and[T^d), (p and ip can be well 
fitted by power law functions of Take (p as an example. 


^^c^ + cPo. 


(105) 


We plot (1 — AT^) in Fig. 18’d) and find that ln(l — AT^) 


can be well fitted linearly with respect to (p, 
ln(l - iF2)«/C + A 


(106) 


Currently we do not have derivations for this linear re¬ 
lation. The good thing is that, in the mass function, 
(1 — K^) is a minor factor. Therefore, as verified in 
Fig.[Tfd) , the mass function can be reduced to 


r 

m = - 
2 

T- 






(107) 




2A 


where a, b, A, and erg are defined in Eqs. (103) and (104). 
We list the fitting results below: 
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Figure 18: (color online). Constraint equati ons and solutions along the slice {x = 0.5, t = t) at the early/slow stage of null 
scattering, (a) and (b): constraint equations (411 and (421. (c) r^t~a{t+by + d, a = (—1.67±0.01) x 10“®, b = —0.6622±0.0008, 
c = 7.897 ± 0.005, d = (-2.610 ± 0.003) x lO’'^. (jRsaln(t + b) + c, a ^ 34.39 ± 0.03, b = -0.245 ± 0.001, c = 3.02 ± 0.04. 
0«a(t + 6)'= + d, a = (-9.7 ± 2.3) x 10'^, b = -0.71 ± 0.03, c = 8.3 ± 0.2, d = -0.133340 ± 0.000005. (d) + bf + d, 

a = (5.0 ± 0.1) X 10■^ b = -0.362 ± 0.006, c = 5.40 ± 0.02, d = (1.19 ± 0.02) x 10“^ ln(l - K^)p 
12.71 ±0.02. Inms 


b = -0.362 ± 0.006, c = 5.40 ± 0.02, d = (1.19 ± 0.02) x 
saln(t ± 6) ± c, a = 77.77 ± 0.02, b = -0.1994 ± 0.0004, c = -15.80 ± 0.03. 


sot ± b, o « —8.73 ± 0.01, 


(i) r t«a(t + by + d, a = (-1.71 ± 0.01) x b = 

(-5.877 ± 0.008) x 10-\ c = 7.880 ± 0.005, d = 
(-2.639 ±0.003) X 

(ii) CT«aln(t ± 6) ± c, a = 34.37 ± 0.03, b = (-1.72 ± 
0.01) X 10-1, c = 3.06 ±0.04. 

(iii) (()«a(t±5)^±d, a = (-2.7±0.6) x b = (-7.6± 
0.3) X 10-1, c = 7.4±0.2, d= (-1.3336±0.0001) x 
10-1. 

(iv) ipKia{t + by + d, a = (5.1±0.1) x 10-^, b = (-2.91± 
0.06)xl0-i, c = 5.38±0.02, d = (1.21±0.02)xl0-3. 

(v) ln(l-7i:2)«at±6, a « -8.73±0.01, b « 12.05±0.02. 

(vi) lnTO«aln(t ± 6) ±c, a = 77.69 ±0.02, b = (—1.268 ± 
0.004) X 10-1, c = -15.65 ± 0.03. 


VII. WEAK SCALAR CHARGE SCATTERING 

In this section, we consider charge scattering with a 
weak scalar field. Parameter settings in this section are 
almost the same as those in the last section with the 
following exceptions: 

(i) Physical scalar field: ^{x,t)\t=o = a ■ 

exp [—(x — Xo)^/b], a = 0.03, 6=1, and Xq = 4. 

(ii) Grid. Spatial range: x G [—12 12]. Grid spacings: 
Ax = At = 0.002. 

As discussed in the above section, in some spacetime 
regions where the scalar field is strong, the inner horizon 
can contract to zero volume, and the central singularity 
becomes spacelike. However, this does not always neces¬ 
sarily happen. After all, it takes energy for the inner 
horizon to contract. When the scalar field carries less 
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Figure 19: Evolutions for charge scattering with a weak scalar field, (a)-(d): evolutions for r, f', and ip. The time interval 
between two consecutive slices is 120At = 0.24. The central singularity is not approached. 


energy, the inner horizon may only contract to a nonzero 
value. This is confirmed by our numerical results plot¬ 


ted in Fig. 19 The evolution of tr is similar to that in 
strong scalar field case and is skipped. These results are 
in agreement with the mathematical proof in Ref. |23j 
and the numerical work in Ref. |34j . Since in this case 
the inner horizon is not totally destructed, one needs to 
reconsider whether the strong cosmic censorship conjec¬ 
ture is valid here. In addition, in Ref. m, it was argued 
that, when Hawking radiation is taken into account, this 
censorship may also be violated. Note that in Ref. m 
the interior of a Schwarzschild black hole was also dis¬ 
cussed with the backreaction from the Hawking radiation 
being taken into account. 

The dynamics for the quantities r, ? 7 , cr, (p, and ip are 
"^e), respectively. The numerical 


plotted in Figs. 20 


results show that, at the late stage, the field equations 
for such quantities become null, in the sense that the 
temporal and spatial derivatives are almost equal, i.e., 
a tt ~ o'^xx- Moreover, the derivatives have oscillations. 
As shown in Fig. [20|[f), the mass function keeps growing 


even as r approaches a constant value. Further details 
are skipped. 


VIII. DARK ENERGY f{R) SINGULARITY 
PROBLEM IN CHARGE SCATTERING 


Up to this point, in the numerical simulations of charge 
scattering that we have implemented in this paper, the 
scalar degree of freedom /' asymptotes to zero as the 
center is approached. Note that, as discussed in Sec. |Vj 
in neutral scattering in dark energy f{R) gravity, when 
the initial velocity or acceleration of /' is large enough, /' 
can go to 1 before the central singularity is approached. 
Consequently, the Ricci scalar R goes to infinity, and the 
simulation breaks down. Next we will show that such a 
problem also happens in charge scattering. 

In the simulation, the values of the parameters are the 
same as those described at the beginning of Sec. |VI[ 
including grid spacings Ax = At = 0.005. How¬ 
ever, the initial value for (p takes the following format: 
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Figure 20: (color online). Dynamics and evolutions on the slice (x = 4, t = t) in weak scattering, (a)-(e): dynamical equations 
for r, rj, a, (j), and ip. (f) evolutions of r and m. 


4>{x,t)\t=o = aexp[-(a: - xo)^] + <^o, with F'(^o) = 0, 
a = —0.05, and Xq = 0.3. We plot the numerical results 
in Fig.|^ As shown in Fig.[2l](b), as the inner horizon is 
approached, /' goes to 1, and the Ricci scalar R becomes 
singular. The simulation breaks down. 


Now we explore the causes of this singularity problem. 
As an example, in Figs. H^d) and [^e), we plot the 
terms in the dynamical equation for (p on the slice 
(x = —2,t = t), and find that initially because of the 
contribution from 4',xx, changes from = 0 to 
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Figure 21: (color online). A singularity problem in charge scattering for the Hu-Sawicki model (30 1 , f{R) = R—DRqR/{R+Ro). 
(a) and (b): evolutions of r and /'. The time interval between two consecutive slices is 20At = 0.01. (c) evolutions of r, /', 
and f't on the the slice (x = —2,t = t). (d) and (e): dynamical equation for (f) on the slice {x = —2,t = t). In (e), as the inner 
horizon is approached, 4>,tt ~ —2r^t4>,t- This equation describes a positive feedback, since —2r,t/r is positive. As a result, when 
is positive, 0 can be accelerated to zero rapidly. Correspondingly, f' goes to 1 as plotted in (b) and (c), and the Ricci scalar 
R becomes singular. Then the simulation breaks down as shown in (a). 
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Figure 22: (color online). Avoidance of the singularity problem in charge scattering for the combined model (1091, f{R) = 
R —DRqR/{R+Ro) + olB?. (a) and (b): evolutions of r and /'. The time interval between two consecutive slices is 60At = 0.15. 
(c) dynamical equation for ij) on the slice {x = —2, t = t). (d) evolution of m on the slice {x = —2, t = t). 


(j)^t > 0 at late time. Near the inner horizon, there 
is 

- ~ - . (108) 

Then cj) is accelerated by gravity. After the inner horizon 
is met, there is ~ (f’.xx > 0. Eventually, 4>, /', and R 
go to 0, 1, and +oo, respectively. The simulation breaks 
down. 

Similar to neutral scattering, the singularity problem 
can be avoided by adding an term to the Hu-Sawicki 
model, 

f{R)=R-§^^+aR\ (109) 

The parameters take the same values as in the last sub¬ 
section with the following exceptions: 

(i) Grid spacings: Aa; = At = 0.0025. 

(ii) Physical scalar field ip: ip{x,t)\t^o = a ■ 
exp [—(a: — Xo^/b ], a = 0.05, 6=1, and xq = 2. 


(hi) Scalar degree of freedom (p: (p{x,t)\t=o = a ■ 

exp[—(a; — a;o)^] + (po, with V'{(po) = 0, a = —0.05, 
and Xq = —2. 


(iv) f{R) model (1091: D = 1.2, Rq = IQ-^, and a = 1. 


The numerical results for charge scattering for this 
modified model are plotted in Fig. As shown in 

Fig. [^b), for this model, f' can cross 1 without dif¬ 
ficulty. The simulation can run smoothly. 


IX. SUMMARY 

In this paper, we studied scalar collapses in flat, 
Schwarzschild, and Reissner-Nordstrdm geometries in 
f{R) gravity numerically. Approximate analytic solu¬ 
tions for different types of collapses were partially ob¬ 
tained. One dark energy f{R) singularity problem was 
discussed. We summarize our work on computational 
and physical issues separately below. 
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A. Computational issues 


(i) The Jordan frame vs the Einstein frame. The field 
equations for f{R) gravity in the Jordan frame 
are more complex than those in general relativity. 
Therefore, for ease of computation, we transform 
f{R) gravity from the Jordan frame into the Ein¬ 
stein frame, in which the formalism can be formally 
treated as Einstein gravity coupled to a scalar field. 

(ii) dudv vs {—dt^ + dx^) in double-null coordinates. In 

the studies of mass inflation, the dudv format of the 
Kruskal-like coordinates, dudvr^diJf, 

is usually used. In the field equations, many terms 
are mixed derivatives of u and v, e.g., In 

this paper, we used the {—dt^ dx^) format in¬ 
stead, ds^ = dt^ -I- dx^) -I- with u = 

{t — x)f 2 = const and v = {t-\- x)/2 = const. In the 
{t,x) line element, one coordinate is timelike, and 
the rest are spacelike. We are used to this setup. It 
is more convenient and more intuitive to use this set 
of coordinates. Moreover, for the {t, x) choice, spa¬ 
tial and temporal derivatives are usually separated, 
e.g., {r^tt - r^xx)- 

We set the initial conditions close to those in a 
Reissner-Nordstrom geometry. With this setup, it 
is convenient to test the code. Removing the terms 
related to the scalar fields, we can test our code 
by comparing the numerical results to the analytic 
ones in a Reissner-Nordstrom geometry. Moreover, 
by comparing numerical results for charge scatter¬ 
ing to the dynamics in the Reissner-Nordstrom geo¬ 
metry, we can obtain intuitions as to how the scalar 
fields affect the geometry. 


(iii) Cauchy horizon: infinite or local regions? As im¬ 
plied by Eq. (28), the exact inner horizon r = t- 
is at the regions where uv and {t^ — x^) are infin¬ 
ite. However, r still can be very close to the inner 
horizon even when uv and (t^ — x^) take moder¬ 
ate values. Consequently, at regions where uv and 
{t^ — x^) take some moderate values, the scalar fields 
and the inner horizon still can have strong interac¬ 
tions, resulting in mass inflation. 


B. Physical issues 

(i) Scalar collapse in f{R) gravity vs scalar collapse 
in general relativity. In scalar collapse, the scalar 
degree of freedom </)(= ^^3/2 In f'/y/SirG) plays a 
similar role as a physical scalar field in general re¬ 
lativity. Regarding the physical scalar field in f{R) 
case, when is negative (positive), the physical 
scalar field is suppressed (magnified) by p. 

(ii) The inner horizon in a Reissner-Nordstrom black 
hole vs the central singularity in a Schwarzschild 
black hole. These two share some similarities. 


Eor Reissner-Nordstrom and Schwarzschild black 
holes, throughout the whole spacetime, the Misner- 
Sharp mass function is constant. When a scalar field 
impacts the inner horizon of a Reissner-Nordstrom 
black hole, the scalar field can modify the geometry 
in the vicinity of the inner horizon significantly, es¬ 
pecially on r t. The inner horizon contracts and 
mass inflation takes place. In neutral scalar col¬ 
lapse toward a Schwarzschild black hole formation, 
the scalar field can also modify the geometry in the 
vicinity of the central singularity dramatically, es¬ 
pecially on the metric component cr |20j . Then mass 
inflation also happens. 

The Belinskii, Khalatnikov, and Lifshitz (BKL) con¬ 
jecture is an important result on dynamics in the 
vicinity of a spacelike singularity [75H75] . The first 
statement of this conjecture is that as the singu¬ 
larity is approached, the dynamical terms dominate 
the spatial terms in the field equations. In other 
words, the way gravity changes over time is more 
important than the variation of the gravitational 
field from one location to the next [79]. We would 
like to say that, to a large extent, later evolutions 
in a strong gravitational field largely erase away 
the initial information on the connections between 
neighboring points. As discussed in Ref. [20] and 
also in this paper, in double-null coordinates, using 
the above argument, one can interpret the follow¬ 
ing behaviors displayed in numerical simulations: 
near the central singularity of a Schwarzschild black 
hole and also near the inner horizon of a Reissner- 
Nordstrom black hole, there are 


i’.x^r^x .. 

V'.t ’’.t 

(IIO) 

>,tt~ - (P,t- 

(III) 


In this paper, it was shown that Eq. (IIO) can ex¬ 
plain the causes of mass inflation, while Eq. (Ill) 
can explain the dark energy /(i?) singularity prob¬ 
lem in collapse. 


The second and third statements of the BKL con¬ 
jecture are that i) the metric terms will dominate 
the matter field terms, while the matter field may 
not be negligible if it is a scalar field; ii) the dy¬ 
namics of the metric components and the matter 
fields is described by the Kasner solution. These two 
statements were confirmed in simulations of neutral 
scalar collapse in f{R) gravity in Ref. [50] and in 
general relativity in Ref. jH] . The second statement 
was also verified in charge scattering in this paper. 
However, the third statement on Kasner solution 
may not apply to the dynamics near the inner hori¬ 
zon in charge scattering. 


(iii) Compact stars vs black holes in f{R) gravity. The 
internal structure of compact stars is usually in an 
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equilibrium state and is static. In f{R) gravity, in¬ 
side compact stars, the scalar degree of freedom, /', 
can be coupled to the energy density of the stars 
and then is not very free to move. However, due 
to strong gravity, the internal structure of black 
holes is dynamical. /' and the matter fields are de¬ 
coupled. As a result, f' is more free to move than 
in the compact stars case. It can keep increasing or 
decreasing until singularities are met. 

(iv) Dark energy f{R) singularity problem: cosmology 
(or static compact objects) vs black hole physics. 
In dark energy f{R) gravity, the Ricci scalar R 
can be singular in both cosmology and black hole 
physics. We consider a homogeneous cosmological 
model. Using the flat Friedmann-Robertson-Walker 
metric, 

ds^ = —dt^ + a^(t)dx^, ( 112 ) 

the equation of motion for f' is 

/" + 3Hf' + U{f) + = 0, (113) 

where H is the Hubble parameter. Due to the finite¬ 
ness of the potential barrier, the force from a per¬ 
turbation of SttT/S may push /' to 1. Correspond¬ 
ingly, the Ricci scalar goes to singularity m In the 
black hole case that we discuss ed i n this paper, the 
equation of motion for 4‘{=y^5/2hif'/\/8 ttG) has 
more complex structure ([M|), 


geometry for the R^ model was simulated. Similar 
results were obtained. Namely, the classical singu¬ 
larity problem, which is present in general relativity, 
remains in collapse in these models. 

(vi) Inside vs outside black holes: local vs global. 
Throughout the whole spacetime of stationary 
Schwarzschild and Reissner-Nordstrom black holes, 
the Misner-Sharp mass function is equal to the black 
hole mass. For a gravitational collapsing system, at 
asymptotic flat regions, the mass function describes 
the total mass of the dynamical system. However, 
in this system, near the central singularity of a 
Schwarzschild black hole or near the inner horizon 
of a Reissner-Nordstrom black hole, the dynamics 
is local. Then the mass function does not provide 
global information on the mass of the collapsing sys¬ 
tem. 

In summary, in this paper, we studied scalar collapses 
in flat, Schwarzschild, and Reissner-Nordstrom geomet¬ 
ries in f{R) gravity. For convenience and intuitiveness, in 
simulating scalar collapses in Schwarzschild and Reissner- 
Nordstrom geometries, Kruskal and Kruskal-like coordin¬ 
ates were used, respectively. Approximate analytic solu¬ 
tions for different types of collapses were partially ob¬ 
tained. Causes and avoidance of a dark energy f{R) 
singularity problem in collapse were discussed. 
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As the central singularity of a Schwarzschild black 
hole or the inner horizon of a Reissner-Nordstrom 
black hole is approached, the above equation can be 
simplified as (111 I, and gravity from the black hole, 
—2r^t4‘,t/’>'■: can cause a similar singularity problem 
as in cosmology or static compact stars. 


(v) The combined f{R) model and the R^ model: sin¬ 
gular or non-singular? At the center of a Schwar¬ 
zschild black hole and in the very early Universe, 
the tidal forces are singular, and general relativity 
fails. Taking into account quantum-gravitational ef¬ 
fects, Starobinsky obtained an model, f{R) = 
R -f aR^. This model has a non-singular de Sitter 
solution, which is unstable both to the past and to 
the future [12 03 ]. In Sec. scalar collapse in a 
Schwarzschild geometry for the combined model (a 
combination of dark energy model and R^ model) 
was explored. A new Schwarzschild black hole, in¬ 
cluding a new central singularity, can be formed. 
Moreover, under certain initial conditions, f' and 
R can be pushed to infinity as the central singular¬ 
ity is approached. In Ref. scalar collapse in flat 
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Appendix: Equations of motion for a physical scalar 
field and f' in the Einstein frame 

In this appendix, based on the equations of motion for 
a massless physical scalar field and the scalar degree 
of freedom f' in the Jordan frame for f{R) gravity, we 
derive the corresponding equations in the Einstein frame. 
In the transformation from the Jordan frame into the 
Einstein frame, we use = X '= \/3/21nx, 
where a tilde denotes that the quantity is in the Einstein 
frame. 

For a scalar field i), the first covariant derivatives of t) 
are equal in two frames, since they are both equal to the 
partial derivative [50] : 

(A.i) 

Then for the first contravariant derivative of ijj, there is 

= X ■ (A.2) 
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For we have 
1 


Combining Eqs. (A.3) and (A.4| yields the equation of 




motion for ij} in the Einstein frame (22). 


1 


VW\ 


■ X 


-4 




• X-" • X • r'^d. 


= x[Oip - g^''df,tljd^{lnx)] 
[2 


(A.3) 


= X I □V' - \/ 3 


In the special case of i/; = x = exp(-\/2/3K(()), there is 


(A.5) 


□x = \/ ^kx^Q<('- 


In the Jordan frame, for a massless scalar field ip, there 
is 


□V’ = 0, 


(A.4) 


Combining Eqs. (§, (1^, (l2§, and ( |A.5[ ) gives the 

equation of motion for (p in the Einstein frame (21). 
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